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1  Introduction 

Our  research  program  was  concerned  with  the  synergy  between  the  variational  partial  differential 
equation  and  the  statistical  approaches  to  problems  in  visual  control.  Moreover,  we  have  given  a 
completely  stochastic  interpretation  of  curvature  driven  flows  for  tracking. 

Vision  is  a  key  sensor  modality  in  both  the  natural  and  man-made  domains.  The  prevalence  of 
biological  vision  in  even  very  simple  organisms,  indicates  its  utility  in  man-made  machines.  More 
practically,  cameras  arc  in  general  rather  simple,  reliable  passive  sensing  devices  which  arc  quite 
inexpensive  per  bit  of  data.  Furthermore,  vision  can  offer  information  at  a  high  rate  with  high  reso¬ 
lution  with  a  wide  held  of  view  and  accuracy  capturing  multi-spectral  information.  Finally  cameras 
can  be  used  in  a  more  active  manner.  Namely,  one  can  include  motorized  lenses  mounted  on  mobile 
platforms  which  can  actively  explore  the  surroundings  and  suitably  adapt  their  sensing  capabilities. 

For  some  time  now,  the  role  of  control  theory  in  vision  has  been  recognized.  In  particular,  the 
branches  of  control  that  deal  with  system  uncertainty,  namely  adaptive  and  robust,  have  been  proposed 
as  essential  tools  in  coming  to  grips  with  the  problems  of  both  biological  and  machine  vision.  These 
problems  all  become  manifest  when  one  attempts  to  use  a  visual  sensor  in  an  uncertain  environment, 
and  to  feed  back  in  some  manner  the  information.  These  issues  were  a  key  thrust  in  our  proposed 
research  program. 

Specifically,  we  have  considered  the  following  problems: 

(i)  Dynamic  Geodesic  Active  Contours  for  Tracking:  Visual  tracking  using  active  contours  is  usu¬ 
ally  accomplished  in  a  static  framework.  The  active  contour  tracks  the  object  of  interest  in  a 
given  frame  of  an  image  sequence,  and  then  a  subsequent  prediction  step  ensures  good  initial 
placement  for  the  next  frame.  This  approach  is  unnatural;  the  curve  evolution  gets  decoupled 
from  the  actual  dynamics  of  the  objects  to  be  tracked.  True  dynamic  approaches  exist,  all  being 
marker  particle  based,  and  thus  prone  to  the  shortcomings  of  such  particle-based  implementa¬ 
tions.  In  particular,  topological  changes  are  not  handled  naturally  in  this  framework.  The  now 
’’classical”  level  set  approach  is  tailored  for  codimension  one  evolutions.  However,  dynamic 
curve  evolution  is  at  least  of  codimension  two  We  have  proposed  a  natural,  efficient  approach 
for  dynamic  curve  evolution  which  removes  the  artificial  separation  of  segmentation  and  pre¬ 
diction,  while  retaining  all  the  desirable  properties  of  level  set  formulations.  This  is  based  on 
a  new  energy  minimization  functional  which  for  the  first  time  puts  dynamics  into  the  geodesic 
active  contour  framework. 

(ii)  Knowledge-Based  Tracking.  One  of  the  key  themes  of  our  work  has  been  the  combination 
of  geometric  partial  differential  equation  methods  with  global  statistics  for  tracking  and  con¬ 
trolled  active  vision.  We  have  developed  a  way  of  combining  statistical  and  anisotropic  dif¬ 
fusion  methods  is  via  a  process  we  call  knowledge-based  segmentation  which  seems  to  give 


some  reasonable  results  for  several  key  applications.  Since  noise  is  in  general  non-additive, 
anisotropic  diffusion  and  related  techniques  directly  applied  to  the  image  do  not  produce  satis¬ 
factory  results.  Moreover,  these  techniques  do  not  introduce  prior  information  about  the  number 
of  objects  present  in  the  scene  when  directly  applied  in  the  image  space.  In  our  approach,  we 
combine  Bayes’  rule  with  either  anisotropic  or  isotropic  diffusion,  introducing  a  priori  knowl¬ 
edge  into  the  segmentation  process  and  solving  the  non-additivity  problem  of  the  noise.  We  also 
extend  this  approach  to  the  segmentation  of  video  data,  incorporating  basic  learning  capabilities 
to  the  knowledge. 

(iii)  Stochastic  Methods  for  Curvature  Driven  Flows :  Many  PDE’s  in  image  processing  and  com¬ 
puter  vision  are  based  on  curvature  driven  flows  from  interfacial  physics.  Such  curvature  flows 
have  been  extensively  considered  from  a  deterministic  point  of  view.  They  have  been  shown 
to  be  useful  for  a  number  of  applications  including  crystal  growth,  flame  propagation,  and 
computer  vision.  In  recent  work,  we  have  described  a  random  particle  system,  evolving  on 
the  discretized  unit  circle,  whose  profile  converges  toward  the  Gauss-Minkowsky  transforma¬ 
tion  of  solutions  of  curvature  driven  flows.  Our  computational  research  program  employs  this 
methodology  as  a  new  way  of  evolving  curves  and  surfaces  for  a  powerful  alternative  to  level 
set  methods.  It  is  a  cornerstone  in  our  computational  program  of  combining  PDE’s  and  statistics 
in  imaging. 

2  Summary  of  Work 

In  this  section,  we  outline  some  key  concepts  from  dynamic  contour  tracking,  and  knowledge-based 
segmentation  methods  upon  which  our  new  research  program  is  based. 

2.1  Dynamic  Contour  Tracking 

Typical  level  set  tracking  techniques  only  make  use  of  static  models  of  active  contours.  In  this  section 
based  upon  our  previous  work  in  [  18],  wc  review  how  dynamics  may  be  put  quite  naturally  into  this 
geometric  framework  giving  us  a  novel  model  of  dynamic  active  contours. 

2.1.1  Parametrized  Dynamic  Curves 

We  consider  the  evolution  of  closed  curves  of  the  form  C  :  S1  x  [0,  r)  »— >  1R2  in  the  plane,  where 
C  =  C(p,  t)  and  C(0,  t)  =  C(l,  t)  [6],  with  t  being  the  time,  and  p  £  [0, 1]  the  parametrization  of  the 
curve.  The  classical  formulation  for  dynamic  curve  evolution  is  derived  by  means  of  minimization  of 
the  action  integral 


(1) 


where  the  subscripts  denote  partial  derivatives.  The  Lagrangian  L  =  T  -  U  is  the  difference  between 
the  kinetic  and  the  potential  energy.  The  potential  energy  of  the  curve  is  given  by 
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where  g  is  some  potential  function  (with  the  desired  location  of  the  curve  forming  a  potential  well); 
UeU  Urig ,  and  Upj  are  the  elasticity,  rigidity  and  potential  field  contributions,  with  their  (possibly 
position-dependent)  scalar  weights  W\,  and  w 2.  A  common  choice  for  the  potential  function  is 


9 ^  1  +  ||G  *  V/(®)||r’  (2) 

where  x  =  [x,y]J  arc  the  image  coordinates,  /  is  the  image,  r  is  a  positive  integer,  and  G  is  a 
Gaussian  of  variance  cr2.  The  kinetic  energy  is 

T  =  J0  ^l|C«|| 2  dp,  (3) 

where  g  corresponds  to  mass  per  unit  length.  The  Lagrangian  used  is  then 

L  =  Jo  \fi\\Ct\\2  -  ^i||Cp||2  -  iHMI2  -  .9(C)  dp.  (4) 

Computing  the  first  variation  SC  of  the  action  integral  (1)  and  setting  it  to  zero  yields  the  Eulcr- 
Lagrange  equations  for  the  candidate  minimizer  in  force  balance  form: 

~  ~~  ~dp ~  ^9' 

This  formulation  is  not  intrinsic  with  respect  to  the  geometry  of  the  curve,  since  it  is  dependent  upon 
its  parametrization,  p. 

2.1.2  Geometric  Dynamic  Curves 
Minimizing  equation  (1)  using  the  Lagrangian 

L  =  Jo  (\lA\ctW2- 9]  \K\\  dp,  (6) 

instead  results  in 

pCtt  =  -p(T  ■  Cts)Ct  -  ,i(Ct  •  Ct,)T  -  (^p\\ct\\2  -  q)kM  -  (V.9  •  (7) 

which  is  intrinsic  and  a  natural  extension  of  the  geodesic  (conformal)  active  contour  approach  [10]. 
Here  J\f  is  the  unit  inward  normal,  and  T  =  the  unit  tangent  vector  to  the  curve,  k.  =  CSH  •  N 
denotes  curvature  and  5  is  arclength. 

Equation  (7)  describes  a  curve  evolution  that  is  only  influenced  by  inertia  terms  and  information 
on  the  curve  itself.  To  increase  robustness,  region-based  terms  could  be  included  in  the  formulation. 
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Normal  Geometric  Dynamic  Curve  Evolution 


To  be  able  to  interpret  the  behavior  of  the  curve  evolution  equation  (7)  it  is  instructive  to  derive  the 
corresponding  evolution  equations  for  the  tangential  and  normal  velocity  components  of  the  curve. 
We  can  write 

Cl  =  or(p,  t)T  +  /3(p,  t)N,  (8) 

where  the  parametrization  p  is  independent  of  time  and  travels  with  its  particle,  and  a  and  (3  corre¬ 
spond  to  the  tangential  and  the  normal  speed  functions  respectively.  By  substituting  equation  (8)  into 
equation  (7),  we  obtain  the  two  coupled  partial  differential  equations: 


at 


Pi 


-(a2)s  +  2na(3, 


—(a{3)s  + 


1  3  o.  1 

y  -2a>  v 


ti - S7q  •  fif. 


(9) 


Clearly,  —  (a2)s  and  —  {a{3)s  are  the  transport  terms  for  the  tangential  and  the  normal  velocity  along 
the  contour,  and  gti  —  V g  •  Jsf  is  the  well  known  geodesic  active  contour  image  influence  term  [10]. 
Note,  that  in  contrast  to  the  static  geodesic  active  contour,  this  term  does  not  directly  influence  the 
curve’s  position,  but  the  curve’s  normal  velocity.  It  resembles  a  force.  Finally,  the  terms  2 na(3  and 
|cv2)k  incorporate  the  dynamic  elasticity  effects  of  the  curve.  If  we  envision  a  rotating  circle 
we  can  interpret  the  term  {\fi2  —  |q;2)ac  as  a  rubberband  (i.e.,  if  we  rotate  the  circle  faster  it  will 
try  to  expand,  but  at  the  same  time  it  will  try  to  contract  due  to  its  then  increasing  normal  velocity; 
oscillations  can  occur).  If  we  restrict  the  movement  of  the  curve  to  its  normal  direction  (i.c.,  if  we  set 
a  =  0)  we  obtain 

A  =  +  ~9K  ~ -Vp  ■  M.  (10) 

2  /i  n 

This  is  a  much  simpler  evolution  equation  In  our  case  it  is  identical  to  the  full  evolution  equation  (9)  if 
the  initial  tangential  velocity  is  zero.  The  image  term  g  only  influences  the  normal  velocity  evolution 
/ 3 .  It  does  not  create  any  additional  tangential  velocity. 

If  there  is  an  initial  tangential  velocity,  and/or  if  the  image  influence  g  contributes  to  the  nor¬ 
mal  velocity  (3  and  to  the  tangential  velocity  a,  the  normal  evolution  equation  will  not  necessarily 
be  equivalent  to  the  full  evolution  equation  (9).  We  can  always  parameterize  a  curve  such  that  the 
tangential  velocity  term  vanishes.  Specifically,  if  we  consider  a  reparametrization 

C(q,  t)  =  C(4>(q,  t),  t),  (11) 


where  <j>  :  R  x  [0,  T)  >  R,  p  =  <j>{q,  t),<j>q  >  0,  then 


0C_dC  dCd$ 
dt  dt  ^  dp  dt 


The  time  evolution  for  C  can  then  be  decomposed  into 

Ct  =  aT  +  f3Af  =  (a(cj)(q,t),t)  +  \\Cp{(f){q,  t),  t)\\<l>t)T  +  /SA/\  (13) 


where 


a  =  a(4>{q,  f),  t)  +  \\Cp(<f>{q,  t),  t)\\</>t, 

P  =  *)>*)■ 
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If  we  choose  0  as 


4>{q,t)t  =  - 


a(0M),O 

IICp^O-OII’ 


(14) 


we  obtain 

Ct  =  (3M,  (15) 

which  is  a  curve  evolution  equation  without  a  tangential  component.  For  all  times,  t ,  the  curve  C 
will  move  along  its  normal  direction  However,  the  tangential  velocity  is  still  present  in  the  update 
equation  for  (3.  After  some  algebraic  manipulations,  we  arrive  at 


li(ftp(pt  -I-  fit)  — 


1 

2^P  +  9 


k.  -  (Vg  Af)Af, 


(16) 


which  depends  on  the  time  derivative  of  the  reparametnzation  function  0,  which  in  turn  depends  on 
the  tangential  component  a.  The  left  hand  side  of  Equation  (16)  represents  a  transport  term  along  the 
curve,  the  speed  of  which  depends  on  the  time  derivative  of  the  reparametrization  function  0. 


2.1.3  Level  Set  Formulation 

There  are  different  ways  to  implement  the  derived  curve  evolution  equations;  see  for  example  [15]. 
We  distinguish  full  and  partial  level  set  implementations.  In  the  full  case,  curves  evolve  in  a  space 
consistent  with  the  dimensionality  of  the  problem.  Geometric  dynamic  curve  evolution  would  thus  be 
performed  in  R1  in  the  simplest  case.  Normal  geometric  dynamic  curve  evolution  would  be  at  least  a 
problem  in  R*.  If  n  is  the  dimensionality  of  the  problem  the  curve  will  be  implicitly  described  by  the 
intersection  of  n  1  hypersurfaces  or  the  zero  level  set  of  an  n-dimcnsional  vector  distance  function. 
Full  level  set  approaches  of  this  form  are  computationally  expensive,  since  it  is  not  obvious  how  to 
devise  a  methodology  comparable  to  a  narrow  band  scheme. 

A  partial  level  set  approach  uses  a  level  set  formulation  for  the  propagation  of  an  implicit  descrip¬ 
tion  of  the  curve  itself  (thus  allowing  for  topological  changes),  but  explicitly  propagates  the  velocity 
information  associated  with  every  point  on  the  contour  by  means  of  possibly  multiple  transport  equa¬ 
tions.  It  trades-off  computational  efficiency  (a  narrow  band  implementation  is  possible  in  this  case) 
for  object  separation:  tracked  objects  that  collide  will  be  merged. 

In  what  follows  we  will  restrict  ourselves  to  a  partial  level  set  implementation  of  the  normal 
geometric  dynamic  curve  evolution.  See  [17]  for  a  detailed  discussion  of  the  full  level  set  method  in 
this  context. 


Partial  Level  Set  Approach  for  the  Normal  Geometric  Curve  Evolution 

The  curve  C  is  represented  as  the  zero  level  set  of  the  function 

V(x(t),t):  R2xR+w|  (17) 

where  x(t)  —  (x(t),  ?/(£))7  is  a  point  in  the  image  plane.  We  assume  '!'>()  outside  the  curve  C  and 
^  <  0  inside  the  curve  C.  Since  the  evolution  of  the  curve’s  shape  is  independent  of  the  tangential 
velocity  we  can  write  the  level  set  evolution  equation  for  an  arbitrary  velocity  Xt  as 

d/t  +  ||Vd'||.A/r  ■  Xt  =  0.  (18) 
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In  our  case  xt  —  pN,  where 


09) 


Equation  (18)  then  becomes 


Wn' 


where  (3  is  given  by  equation  (10).  Using  the  relation 


k  =  V- 


(20) 

(21) 


we  obtain 


Pt  —  ( -0 2  H — q)k  H — Vg  • 

2  n  p 


Wn' 


(22) 


We  need  to  propagate  the  normal  velocity  P  along  with  the  zero  level  set  of  \P.  This  could  be  accom¬ 
plished  by  the  transport  equation 


Pt=P 


1 Wn 


■vp, 


(23) 


which  is  not  natural,  since  equations  (22)  and  (23)  are  both  first  order  evolution  equations  for  f3. 
Instead  of  using  equation  (23)  we  can: 


(a)  Neglect  equation  (23),  but  compute  extension  velocities  at  every  time  step.  Sinee  the  extensions 
are  normal  to  the  contours,  normal  propagation  of  the  level  set  function  will  guarantee  a  constant 
velocity  value  along  the  propagation  direction  (up  to  numerical  errors).  Specifically  V/3±V^ 
in  this  case  and  thus 

•  V/3  =  0.  (24) 


(b)  Incorporate  equation  (23)  into  equation  (22). 

To  accomplish  the  latter  we  change  our  Lagrangian,  and  extend  it  over  a  range  of  level  sets.  For  each 
time  t ,  and  0  <  r  <  1  let 

c(r)(f)  :=  {(x,y)  E  R2  =  r}.  (25) 


Using  the  Lagrangian 

L- 111.0**-°) 

|  ds  dr, 

(26) 

we  obtain  the  action  integral 

‘-III.  (>■-•) 

|  ds  dr  dt, 

(27) 

which  is 

C 

=  [  [  [ (\vP2  -  9 )  ds  dt  dr 

Jo  Jo  Jc  1 

=  L  (.Jo  JA^  ~~  9 ^  ^  ^  dr )  dt 

=  /  ( {\nP2  ~  9)\\W\\dx  dy  dt, 

Jo  J  n  2 

(28) 
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where  TL 1  is  the  one-dimensional  Hausdorff  measure  and  we  applied  the  coarea  formula.  This  casts 
the  minimization  problem  into  minimization  over  an  interval  of  level  sets  in  a  fixed  coordinate  frame 
( x  and  y  are  time  independent  coordinates  in  the  image  plane).  Using  equation  (20)  we  can  express  0 
as 


(29) 


Substituting  (29)  into  equation  (28)  yields 


C  = 


2||V*|| 


<?||V\I/||)  dx  dy  dt  :=  S  4'], 


(30) 


which  is  the  new  ^-dependent  action  integral  to  be  minimized.  Then,  6S  —  0  if  and  only  if 


a  \liv*||/  VVf  liv*|iv  ||V*||J 

The  curve  evolution  is  thus  governed  by  the  system: 

(it 

Expanding  equation  (32)  yields 


=  V  •  [  - (-  +  -0 2) 


1  o  1  1  V4> 

A-(^  +  ?)«  +  -V,w  +  Ww, 


(31) 


(32) 


(33) 


which,  as  desired,  combines  the  velocity  update  equation  (22)  with  the  velocity  propagation  equa¬ 
tion  (23).  Naturally,  this  equation  reduces  to  the  result  obtained  from  our  original  minimization 
problem  (see  equation  (22))  if  extension  velocities  are  used  to  tie  the  evolution  of  all  level  sets  to 
the  evolution  of  the  zero  level  set.  Note,  that  the  system  of  equations  (32)  constitutes  a  hyperbolic 
conservation  law  for  the  velocity  v.  The  propagation  of  the  level  set  function  d/  is  described  (as  usual) 
by  an  equation  of  Hamilton-Jacobi  type. 


2.2  Statistics  and  Curvature  Driven  Flows:  Knowledge-Based  Segmentation 

We  summarize  here  some  of  our  results  on  the  synergy  of  geometric  partial  differential  equation  and 
statistical  methods  based  on  knowledge-based  segmentation.  We  combine  Bayes’  rule  with  either 
anisotropic  or  isotropic  diffusion,  introducing  a  priori  knowledge  into  the  segmentation  process  and 
solving  the  non-additivity  problem  of  the  noise.  We  also  extend  this  approach  to  the  segmentation  of 
video  data,  incorporating  basic  learning  capabilities  to  the  knowledge. 

This  has  a  number  of  applications  including  SAR  imagery  [7]  and  missile  tracking  [9,  23].  In 
particular,  in  the  case  of  tracking  missiles  for  the  airborne  laser  (ABL)  program,  wc  are  interested  in 
tracking  the  location  of  the  nose  of  the  missile.  Thus  we  only  need  to  separate  the  relevant  portion  of 
the  missile  from  the  background.  Because  of  the  noisy  nature  of  the  images  due  to  atmospheric  effects 
simpler  thresholding  techniques  (e.g.,  histogramming  the  pixel  distributions  and  trying  to  separate 
the  peaks)  do  not  work  very  well  for  the  missile  videos.  On  the  other  hand,  we  have  found  that 
the  knowledge-based  approach  of  [9,  23]  gave  excellent  results  in  comparison  to  weighted  centroid 
techniques  especially  when  combined  with  geometric  active  contours. 
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2.2.1  Basic  Set-up 


Our  set-up  begins  with  the  assumption  that  the  image  is  composed  of  n  classes  of  objects.  For  se¬ 
quences  of  images,  this  value  n  is  assumed  constant  For  the  missile  problem  [23]  ,  we  assume  two 
classes,  corresponding  to  the  missile  and  the  background.  Thus  in  this  setting  we  will  see  that  wc 
have  a  form  of  adaptive  thresholding.  This  is  the  case  wc  will  describe  here.  The  technique  however 
is  general  and  can  be  applied  to  any  number  of  classes.  (In  [7]  for  SAR  data,  the  number  of  classes 
was  three.)  The  goal  of  our  segmentation  is  to  determine  to  which  class  each  pixel  in  each  image 
belongs.  Our  basic  model  assumes  that  the  value  of  each  pixel  in  a  given  class  can  be  thought  of  as 
a  random  variable  with  a  known  distribution,  and  that  these  variables  are  independent  across  pixels 
(this  just  simplifies  the  exposition,  but  any  distribution  can  be  used).  Thus,  for  example  for  the  case 
of  normal  distributions  the  likelihood  of  a  particular  pixel  i  having  a  certain  value  v  given  that  it  is  in 
class  c  €  {missile,  background}  is: 

Pr(V<  =  v\Ci  =  c)  =  2.  exp  (~\— — j^~)  ,  (34) 

V2nac  \  2  aj  ) 

where  i  is  an  index  ranging  over  all  pixels  in  the  image,  V{  is  the  value  of  the  pixel,  and  Cx  is  its  class. 
As  usual,  //.c  and  rrc  denote  the  mean  and  standard  deviation  of  class  c.  In  practice,  these  parameters 
are  estimated  from  a  set  of  sample  images  or  learned  from  past  frames  for  video  data. 

Next,  we  assume  that  there  is  some  known  prior  probability  that  a  particular  pixel  will  belong  to 
a  certain  class.  For  single-image  data  sets,  we  assume  a  homogeneous  prior,  i.e.,  that  Pr(Cj  =  c)  is 
the  same  over  all  spatial  indices  i  It  is,  however,  possible  to  incorporate  a  priori  knowledge  about  the 
image  here,  for  example  if  it  were  known  that  the  missile  is  more  likely  to  be  near  the  center  of  the 
image  than  near  the  edge.  For  sequences  of  images,  we  have  used  a  learned  prior,  as  described  below. 

Given  a  set  of  intensity  distributions  Pr(Vri  =  v\ Cx  —  c)  and  priors  Pr(C{  =  c),  we  can  apply 
Bayes’  Rule  from  elementary  probability  theory  to  calculate  the  posterior  probability  that  a  given 
pixel  belongs  to  a  particular  class,  given  its  intensity: 


Pr  (Q  =  c\Vi  =  v)  = 


Pr (VS  =  v\Ci  =  c)  Pr(Ci  =  c) 
Pr(Vi  =  v\Ci  =  7)  Pr(Cf  =  7) 


(35) 


Our  approach  is  to  calculate  the  posteriors  Ff  \=  Pr(Ci  —  c\Vt  —  v)  using  the  given  distributions 
and  (35)  above,  and  then  to  apply  either  isotropic  or  anisotropic  smoothing  to  each  Pc  (note  that  the 
denominator  is  just  a  normalization  constant  that  can  be  “ignored”).  Specifically,  we  have  chosen  to 
smooth  by  evolving  Pc  according  to  a  discretized  version  of  the  partial  differential  equation 

pipe 

—  =  {(PCy?PCxx  -  2 PCxPCyny  +  (PtfFyy)1'3-  (36) 

This  equation  defines  the  affine  geometric  heat  flow,  under  which  the  level  sets  of  Pc  undergo  affine 
curve  shortening.  This  particular  diffusion  equation  was  chosen  because  of  its  affine  invariance,  be¬ 
cause  it  preserves  edges  well,  and  because  of  its  numerical  stability  and  ease  of  computation.  The 
goal  of  this  process  is  to  diffuse  information  from  one  pixel  to  the  other,  making  then  a  region-based 
(adaptive)  decision  and  not  just  a  local  pixel-wise  one. 

The  segmentation  is  then  obtained  using  the  maximum  a  posteriori  probability  estimate  after 
anisotropic  smoothing.  That  is. 


c: 


ce{ 


arg  max 

missile,  background 


Pr*(Ci  =  c\Vi  =  v) 


where  Pr*  (Ci  =  c|  =  v)  is  the  smoothed  posterior  probability. 


(37) 
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2.2.2  Application  to  Visual  Tracking 

When  segmenting  sequences  of  images,  we  have  extended  the  model  so  that  information  from  one 
frame  is  used  in  the  segmentation  of  the  next,  and  in  this  way  have  introduced  a  kind  of  learning 
into  the  method  There  are  a  number  of  ways  in  which  this  can  be  done.  By  far  the  most  effective 
way  we  have  found  is  to  modify  our  assumption  of  homogeneous  priors.  In  particular,  we  have 
used  the  smoothed  posteriors  Pc  from  one  frame  as  priors  Pr(Cf  =  c)  in  the  segmentation  of  the 
next  frame.  Thus  this  scheme  fits  into  a  Bayesian  tracking  framework.  We  have  also  tested  relaxing 
our  assumption  that  the  pixel  intensities  are  distributed  according  to  fixed  normal  distributions.  We 
estimated  the  distribution  parameters  of  the  normal  distributions  from  frame  to  frame  by  calculating 
new  sample  means  and  variances  based  on  the  segmentation  of  earlier  images.  Finally,  we  completely 
removed  the  assumption  that  the  intensities  are  normally  distributed  This  was  done  by  calculating 
the  sample  distribution  of  intensities  within  each  class  as  images  were  segmented,  and  then  using  this 
distribution  as  Pr(V*  =  v\C{  =  c )  in  (34)  when  segmenting  succeeding  frames. 

This  type  of  approach  combined  w  ith  active  contours  was  found  to  be  extremely  effective  for  ABL 
missile  tracking  [23].  For  other  scenarios,  it  is  possible  to  introduce  multi-scale  texture  models  for  the 
likelihood  of  the  background.  Another  possible  extension  would  be  to  consider  that  n,  the  number  of 
classes  in  the  image,  is  not  given  and  needs  to  be  estimated  as  well.  This  can  be  done  for  example  via 
expectation  maximization  (EM)  type  algorithms.  Note  though  that  since  the  scheme  here  described  is 
extremely  fast,  especially  for  video  data  where  the  number  of  smoothing  steps  is  dramatically  reduced, 
a  brute-force  search  for  n  in  a  given  range  might  be  good  enough  for  a  number  of  applications. 

This  knowledge-based  approach  has  a  very  strong  connection  to  Markov  random  field  (MRF) 
techniques.  Indeed,  the  prior  smoothing  of  the  posterior  probabilities  gives  the  MAP  solution  to  a 
discrete  MRF  with  a  non-interacting,  analog  discontinuity  field.  The  combination  of  a  discontinuity 
field  with  a  discrete  MRF  can  have  some  important  consequences  since  it  allows  the  disabling  of 
clique  potentials  across  discontinuities.  This  is  in  contrast  to  the  isotropic  (linear)  smoothing  of  the 
posterior  probabilities,  which  corresponds  to  computing  the  MAP  solution  of  a  single  discrete  MRF 
using  continuous  relaxation  labeling.  Thus  we  see  that  the  random  field,  anisotropic  diffusion,  and 
curvature  driven  flow  approaches  to  segmentation  and  tracking  are  intimately  related. 

2.3  Interacting  Particle  Systems  and  Approximate  Flows 

In  typical  methods  such  as  those  described  above  for  knowledge-based  segmentation  for  combining 
statistics  with  PDE’s,  the  underlying  PDE  model  remains  deterministic  and  is  implemented  as  such, 
for  example  using  level  sets.  Wc  have  developed  a  very  different  class  of  algorithms  in  which  the 
PDE  is  eliminated  and  the  evolution  takes  place  using  a  random  process.  We  should  note  that  this 
methodology  first  proposed  for  curvature  driven  flows  in  [2,  3]  is  mathematically  well  founded,  and 
may  be  justified  using  the  theory  of  hydrodynamical  limits. 

The  new  class  of  algorithms  is  very  much  connected  to  scale  in  physical  systems.  The  typical  way 
of  modelling  a  physical  system  is  on  a  macroscopic  scale,  that  is  via  a  continuous  partial  differential 
equation.  Level  set  computational  methods  fall  into  this  category  of  algorithms.  Here  wc  arc  looking 
at  a  microscopic  approach:  collections  of  discrete  particles  undergoing  random  walks.  What  wc  arc 
doing  then  is  replacing  the  PDE  by  an  underlying  microscopic  interacting  particle  system. 

The  theory  of  interacting  particle  systems  goes  back  at  least  to  [22]  who  studied  new  types  of 
random  walk  models  with  interactions  among  the  particles.  There  are  many  examples  including  ex¬ 
clusions  processes,  branch  processes,  and  contact  processes.  There  arc  several  advantages  to  this  type 
of  methodology  including: 
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1  PDE’s  and  continuous  analysis  become  unnecessary. 

2.  No  discretizations  arc  necessary.  The  interacting  particle  system  is  already  constructed  on  a 
discrete  fixed  lattice. 

3.  Increased  robustness  to  noise,  and  ability  to  include  processes  into  the  given  system. 

In  summary,  this  framework  gives  a  completely  stochastic  way  of  modelling  the  key  flows  used  in 
controlled  active  vision  will  all  the  possible  benefits  just  described. 

2.3.1  Basic  Definitions  of  a  Particle  System 

In  the  theory  of  interactive  particle  systems,  particles  are  distributed  on  a  d-dimensional  discrete 
integer  lattice  .  The  nodes  of  the  lattice  represent  sites  where  particles  reside.  The  state  of  the  lattice 
is  the  distribution  of  particles  over  the  sites.  This  state  can  evolve  over  time  due  to  events  such  as 
births  and  deaths  of  particles,  as  well  as  jumps  of  particles  to  neighboring  sites.  The  probability  of  an 
event  happening  at  a  particular  site  can  be  influenced  by  the  number  of  particles  at  neighboring  sites, 
therefore  adding  interaction  between  particles.  If  the  rate  of  events  at  a  site  is  only  influenced  by  the 
number  of  particles  at  that  site,  then  the  process  is  called  a  zero  range  process.  For  our  purposes,  we 
are  interested  in  zero  range  processes  where  particles  only  interact  with  particles  sitting  on  the  same 
site. 

2.3.2  Geometric  Interpretation  of  a  Particle  System 

In  our  AFOSR  sponsored  research,  we  were  able  to  construct  an  interacting  particle  system  that 
evolves  microscopically  on  a  ID  lattice  according  to  the  macroscopic  behavior  of  the  mean  curva¬ 
ture  flow  PDE.  This  is  a  fundamental  result  that  links  the  microscopic  scale  (particles  evolving  inde¬ 
pendently  of  particles  at  neighboring  sites)  to  the  macroscopic  scale  (curve  evolving  according  to  a 
macroscopic  PDE). 

Before  we  detail  our  particle  system,  we  explain  its  geometric  interpretation.  This  interpretation  is 
fundamental  for  the  use  of  particle  systems  as  a  tool  for  stochastic  curve  evolution.  What  we  want  to 
do  is  to  be  able  to  represent  the  particle  system  graphically  as  a  curve.  As  the  particle  system  evolves, 
the  corresponding  curve  is  updated,  resulting  in  a  stochastic  curve  evolution. 

Our  particle  system  is  a  ID  line  of  sites  with  connected  ends,  that  we  call  a  discrete  torus  where 
each  site  is  a  point  space  evenly  on  the  unit  circle.  Therefore  each  site  is  parameterized  by  an  angle. 
Each  site  contains  a  certain  positive  number  of  particles,  denoted  with  the  positive  function  m(6) 
where  6  parameterizes  each  site.  In  order  to  map  our  particle  system  to  a  curve,  we  can  connect 
segments  oriented  according  to  the  angle  of  each  site  and  scaled  by  the  number  of  particles  at  that  site. 

2.3.3  Birth  and  Death  Zero  Range  Particle  Systems 

We  first  set-up  some  notation. Let 


C(p,  t)  :  Sl  x  [0.  T]  —  R2 


be  a  family  of  embedded  curves  where  p  parametrizes  the  curve  and  t  the  family.  T  hen  we  consider 
curvature-driven  flows  of  the  form 


—  =  v(K(p,t))Ar, 
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(38) 


where  k  denotes  the  curvature  and  J\f  the  inward  unit  normal.  Using  the  standard  angle  parametriza- 
tion  0 ,  we  interpret  m(9 ,  t)  :=  1  /k(0,  t)  as  a  density,  and  compute  it  evolution  to  be: 


dm(0 ,  £) 
^~dt 


d2V(m(0,  t)) 

dt ? 


U(x)  r=  V^l/x). 


(39) 


Note  that  V^x)  =  1/x  corresponds  to  the  Euclidean  curve  shortening  flow.  In  Equation  (39),  the 
first  term  on  the  right  hand  side  is  called  the  diffusion  term  and  the  second  term  the  reaction  term. 

Our  interest  is  in  constructing  stochastic  approximations  to  the  solutions  of  the  Equation  (38)  or 
equivalently  Equation  (39).  Approximations  corresponding  to  polygonal  curves  have  been  discussed 
in  the  literature  under  the  name  “crystalline  motion.”  Our  approach  is  different  and  can  be  thought  of 
as  a  stochastic  crystalline  algorithm:  we  construct  a  stochastic  particle  system  whose  profile  defines  an 
atomic  measure  on  S1  ,  such  that  the  corresponding  curve  is  a  polygon,  Applying  standard  tools  from 
hydrodynamic  limits,  it  is  proven  in  [2]  that  the  (random)  evolution  of  this  polygonal  curve  converges, 
in  the  limit  of  a  large  number  of  particles,  to  curve  evolution  under  the  curve  shortening  flow. 

The  approximations  we  use  are  based  on  so-called  birth  and  death  zero  range  particle  systems. 
To  get  a  flavor  of  the  simplicity  of  the  algorithm,  we  write  down  this  system  down  in  some  detail.  Let 
Tj\[  =  Jj/NJj  denote  the  discrete  torus.  Let  g  :  N  — »  R+  (the  jump  rate ,  with  <?(0)  =  0),  6  :  N  — *  R + 
(the  birth  rate),  d  :  N  — *  R+  (the  death  rate ,  with  d( 0)  =  0)  be  given,  and  define  the  Markov 
generator  on  the  particle  configuration  E ^  =  NT*  by 

(CNf)(v)  =  N2(C0f)(ri)  +  (A /)(»?) ,  /  €  Cb(EN) , 


where 


and 


(Co /)(!?)  =  \  £  9(v(i))  [f(ti+1)  +  fW'*-1)  -  2 /(»?)] 

(CifM  = 

£  [&(»?(*))  [f(vl,+  )  -  /(»?)]  +  d(i](i))  [f(rf~)  -  f(rj)]]  , 

f  VU)  +  1,  j  =  i±  1  ,r)(i)  ^  0, 
n,  '±1U)={  ’iti)  - i,  j  =  o, 

[  else 


rru) 


(  ij(j)  +  1  ,j  =  i, 
\  i](j),e Ise 


vl-(i) 


(  vU)  - 1 ,3  =  i,n(i)  >  0, 
\  vti),  else 


Note  that  the  zero-range  part  Co  approximates  diffusion  term  of  equation  (39)  while  the  birth-death 
part  C\  approximates  the  reaction  term  of  (39). 
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2.3.4  Example:  Approximate  Euclidean  curvature  flows 

As  an  example,  in  this  section  we  give  the  stochastic  approximation  of  the  Euclidean  curvature  flow. 
For  the  general  curve  shortening  case;  see  [2,  3]  for  the  details.  We  now  present  candidates  for 
the  functions  b,d,g  defining  the  particle  systems.  Indeed,  for  curve  shortening  (V(x)  =  l/x  in 
Equation  (39)  above)  one  may  show  that  the  rates  are 

g( l)  =  e-2,  g(k)  =  e~1k/(k  -  1),  k  >  2, 

6(2)  =  2e-2,  b(k)  =  0,k^  2, 

d(  l)  =  e-2,  d(k)-0,k^l.  (40) 

This  means  that  as  the  number  of  cites  N  goes  to  infinity  and  as  6  goes  to  zero,  the  stochastic 
particle  system  converges  to  the  solution  of  curve  shortening. 

2.3.5  Stochastic  Snakes 

Wc  have  formulated  a  geometric  active  contour  model  stochastically,  and  derived  a  complete  stochas¬ 
tic  snake  formulation.  In  this  research  program,  we  plan  to  derive  fully  stochastic  versions  of  all 
the  segmentation  results  described  above.  For  the  geometric  active  contour  model  case  which  we 
considered  in  this  preliminary  study,  the  density  function  evolves  according  to 

mt  =  -(m~')ee<t>  -  2 (m_1)e0e  -  m_1(006,  +  0),  (41) 

where  0  is  the  (conformal)  stopping  term,  and  subscripts  indicate  partial  derivatives.  The  correspond¬ 
ing  birth-death  zero  range  process  is  defined  exactly  as  above  with  the  addition  of  the  following  for 
the  first-order  drift  term  on  the  right  hand  side  of  (41): 

N  2(T7(i))(/(r?I,t+1)  -  /(»?))• 
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